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Abstract 

A new method is presented for solving Poisson's equation inside an open- 
ended rectangular pipe. The method uses Fast Fourier Transforms (FFTs) 
to perform mixed convolutions and correlations of the charge density with 
the Green function. Descriptions are provided for algorithms based on the 
ordinary Green function and for an integrated Green function (IGF). Due 
to its similarity to the widely used Hockney algorithm for solving Poisson's 
equation in free space, this capability can be easily implemented in many 
existing particle-in-cell beam dynamics codes. 

Keywords: convolution, correlation, FFT, Poisson equation, Green's 
function, Hockney method 



1. Introduction 

The solution of Poisson's equation is an essential component of any self- 
consistent beam dynamics code that models the transport of intense charged 
particle beams in accelerators, as well as other plasma particle-in-cell (PIC) 
codes. If the bunch is small compared to the transverse size of the beam 
pipe, the conducting walls are usually neglected. In such cases the Hockney 
method may be employed [U EJ |3] . In that method, rather than computing 
N£ point-to-point interactions (where N p is the number of macroparticles), 
the potential is instead calculated on a grid of size (2N) d , where N is the 
number of grid points in each dimension of the physical mesh containing the 
charge, and where d is the dimension of the problem. Using the Hockney 
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method, the calculation is performed using Fast Fourier Transform (FFT) 
techniques, with the computational effort scaling as (2N) d (log22N) d . 

When the beam bunch fills a substantial portion of the beam pipe trans- 
versely, or when the bunch length is long compared with the pipe transverse 
size, the conducting boundaries cannot be ignored. Poisson solvers have been 
developed previously to treat a bunch of charge in an open-ended pipe with 
various geometries jU [5]. Another approach (employed, e.g., in the Warp 
code |6j), is to use a Poisson solver with periodic, Dirichlet, or Neumann 
boundary conditions on the pipe ends, and to extend the pipe in the simula- 
tion to be long enough so that the field is essentially zero there. 

Here a new algorithm is presented for the open-ended rectangular pipe. 
The new algorithm is useful for a number of reasons. First, since its structure 
is very similar to the FFT-based free space method, it is straightforward 
to add this capability to any beam dynamics code that already contains 
the free space solver. Second, since it is Green-function based, the method 
does not require modeling the entire transverse pipe cross section, i.e., if the 
beam was of small transverse extent one could instead model only a small 
transverse region around the axis. Third, since it is based on convolutions and 
correlations involving Green functions, the method can use integrated Green 
function (IGF) techniques. These techniques have the potential for higher 
efficiency and/or accuracy than non-IGF methods [7J E], and are used in 
several beam dynamics codes including IMPACT, MaryLie/IMPACT, OPAL, 
ASTRA, and BeamBeam3D P \M EH 021 E3] . 

The solution of the Poisson equation, V 2 = — p/eo, for the scalar poten- 
tial, </>, due to a charge density, p, can be expressed as, 

(f)(x,y,z) = J J J dx'dy'dz'p(x',y',z')G(x,x',y,y',z,z'), (1) 

where G(x,x',y,y',z,z') is the Green function, subject to the appropriate 
boundary conditions, describing the contribution of a source charge at loca- 
tion (x',y',z') to the potential at an observation location (x,y,z). For an 
isolated distribution of charge this reduces to 




<P{x,y,z)=l I I dx'dy'dz'p(x',y',z')G(x -x',y -y',z- z'), (2) 

where 
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G(u, v, w) = —=====. (3) 
Vii + v z + w z 

A simple discretization of Eq. ([2j on a Cartesian grid with cell size (h x , h y , h z ) 
leads to, 



— h x h y h z ji fc'Gi-ii t k-k' , (4) 

i'=i j'=i k'=i 

where pij^ and Gi-i/j-j/^-k' denote the values of the charge density and the 
Green function, respectively, defined on the grid. 

As is well known [T3], FFT's can be used to compute convolutions by ap- 
propriate zero-padding of the sequences. A proof of this is shown in the Ap- 
pendix, along with an explanation of the requirements for the zero-padding. 
As a result, the solution of Eq. (Ej) is then given by 

o,.jm = h x h y h z F bbb {{pff Phhk ){pffG^ k )} (5) 

where the notation has been introduced that denotes a forward FFT 

in all 3 dimensions, and J rbbb denotes a backward FFT in all 3 dimensions. 
The treatment of the open rectangular pipe relies on the fact, as described 
in the Appendix, that the FFT-based approach works for correlations as well 
as for convolutions, the only difference being the direction of the FFTs. As 
an example, consider for which the y variable involves a correlation 

instead of a convolution. Then, 



7 1 k 

'"max Jmax ^max 

/v.,'./,<^ ;<„././, /,- = J r ' ,Ih {iP ff i>,,jMP' ,f G, J M)}. (6) 

i'=i j'=i k'=i 

Because of this, and the fact the the Green function for a point charge in an 
open rectangular pipe is a function of (x ± x', y ± y', z ± z'), an FFT-based 
algorithm follows immediately. 



2. Poisson's equation in an open rectangular pipe 

The Green function for a point charge in an open rectangular pipe with 
transverse size (0, a) x (0, b) can be found using eigenfunction expansions 
(see, e.g., [IS]). It can be expressed as, 
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. . . 1 ^-^ ^-^ 1 . mux . mux . niry . mry , ,, 

G{x, x , y, y , z, z ) = - — - > > sm sm sm — — sm — — e mn| '. 

m=l n=l mn 

(7) 

where = ( r ^ E ) 2 + (^f 1 ) 2 [16J. Equivalently, it can be expressed as a 
function of a single rectangular pipe Green function, 

G = R(x — x',y — y', z — z') — R(x — x',y + y', z — z') 

—R(x +x',y — y', z — z) + R(x + x',y + y', z — z'), (8) 

where 

R(u,v,w) = -> > cos cos——e KmnM . (9) 

m=l n=l 

It follows that an algorithm for solving the Poisson equation in an open 
rectangular pipe is given by 

4>i,j,k/{h x h y h z ) = 

-^Vftj,*)^'^)} + P Sb {{P Sf p l , 3 , k ){T bbf R l , 3 , k )} (10) 

At each step of a simulation charge is deposited on a doubled grid, and 
its forward FFT, p%j,k), is computed. The function R is tabulated in 

the doubled domain, and 4 mixed Fourier transforms are computed, namely 
(F f,f Rij,k), (J r6// J Rij, fe ), and (J&fRidk)- The transformed 

charge density is multiplied by each of the transformed Green functions, 4 
final FFTs are performed, and the results are added to obtain the potential. 

Note that, although the Green function in Eq. ^ contains 4 terms, only 
a single function R needs to be tabulated on a grid, and it is this tabulated 
function that is transformed in 4 different ways and stored. 



3. Integrated Green Function 

The discrete convolution, Eq. Q, is a simple approximation to Eq. (J2|. 
It makes use of the Green function only at the grid points, even though it is 
known everywhere in space. This can lead to serious inaccuracy when p and 
G have a disparate spatial variation, especially when the grid cells have high 
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aspect ratio. This can happen in a variety of situations (isolated, open-ended 
pipe, etc.) depending on the problem geometry and the number of grid cells 
in each dimension. For example, a simulation of a disk-shaped galaxy would 
have highly elongated cells unless many more cells were used in the wide 
direction. Also, inside a conducting pipe, the beam might be long and slowly 
varying, whereas the pipe Green function decays exponentially with z. 

Qiang described a method for approximating Eq. ^ to arbitrary ac- 
curacy using the Newton-Cotes formula [17) . Integrated Green functions 
(IGF's) provide another means to approximate Eq. ^ accurately when cer- 
tain integrals involving the Green function can be obtained analytically [TJ [S] . 
This is accomplished by assuming a simple analytical form for the variation of 
p within a cell, and analytically performing the convolution (or correlation) 
integrals for each cell of the problem. As a result, the accuracy is controlled 
by how well the discretization resolves p, not G. 

For illustration, consider the 2D problem with isolated boundary condi- 
tions. If p were assumed to be a sum of delta functions at the grid points, 
we would obtain the 2D analog of Eq. Q. Instead suppose that linear basis 
functions are used to approximate p within each cell. Then 

<f>(xi,yj) = 

Pi',3 1 J J dxdy — G{Xi - x v - x , y 3 - y 3 
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"x>"y 



sr f hx f hy a >a ' x ^hLzyl r< ' '\ 

+ y. Pi'+ij' / I dxdy — — G{Xi - x v - x , y 3 - - y r - y ) 

v ., Jo Jo n x n y 

ST f kx f kv A 'A ' (h^K n( > * 

+ 2^ P i 'j'+ 1 I I dxdy — — G{xi - x v - x , y 3 - - y r - y ) 

i>ji Jo Jo n x n y 

Ef hx f hy x'y' 
Pi'+i,f+i / / dx'dy' -r—r-G(xi - x v - x' ,y 3 - - y f - y'). (11) 
j/j, Jo Jo n xii y 

Shifting the indices in the last 3 sums above, collecting terms, and factoring 
out h x h y (in analogy with Eq. Q), the IGF is the coefficient of h x h y pi'ji , 



(f>(xi, yj) = h x h y Pi',j>Gt-i>,j- 3 >- (12) 



1 ->3 



To simplify the discussion that follows, suppose that p is zero on the bound- 
ary; this allows one to shift indices and collect terms without worrying about 
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possible boundary issues. 

Returning to the problem of the open-ended rectangular pipe, here we 
will use the IGF approach to treat the longitudinal dimension. Suppose 
that, inside the k th cell, the longitudinal dependence of the charge density is 
given by, 

p(z) = y \Pk{h z ~(z- z k )) + p k+1 (z - z k )\ . (13) 
The longitudinal dependence of the IGF is therefore, 

9* = Ti I dz ' fafrz - ( z> - z k)) + Pk+i {z - z k )\ e-"™\*- z '\. (14) 

Performing the integrations and summing over all the cells, two terms from 
adjacent cells contribute to the the IGF, with the result, 



1 



9z{w) 



z mn 



2h z K mn S W)0 + (e-^-l^+^l - 2e~ KmnH + e -"™>- fc »l) 



(15) 

where w = z — z k . In summary, in analogy to Eq. ([9]), the integrated Green 
function, Ri nt , integrated in just the longitudinal coordinate, for a distribu- 
tion of charge in an open-ended rectangular pipe is given by, 



oo oo 



1 ^-^ 1 mini nnv 
R int {u, v, w) = - — - > > cos cos ——g z {w). (16) 

m=l n=l 



4. Numerical Example 

Consider a rectangular waveguide of full width and height a = b = 4 cm. 
The Green functions, Eq. ^ and Eq. (16), will be calculated using 



m, n 



1, . . . , 20. Consider a 3D Gaussian charge distribution with transverse rms 
sizes a x = 0.15a = 6mm, a y = 0.156 = 6mm, and longitudinal rms size 
a z . Three cases with different rms bunch length will be considered, a z = 
1.2 cm, a z = 12 cm, and a z = 1.2 m. The distribution is set to zero at 
x 2 /al + y 2 /a 2 + z 2 /a 2 > 3 2 . Fig. [TJ left, shows the charge density and 
Green function as a function of z, down the center of the pipe, for the case 
<j z = 1.2 cm. As shown in the Appendix, the convolution can be performed 
using FFTs by zero-padding the charge density over the domain of the Green 
function, and by treating the Green function as a periodic function. 
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Fig. [TJ right, shows how the charge density and Green function are actu- 
ally stored in memory, using 256 longitudinal grid points: p is left in place 
and zero-padded, and G is circular-shifted so that G(x = 0, y = 0, z = 0) is 
at the lower bound of the Green function array For the case a z = 12 cm, 
Fig. [TJ right, would look similar except that the Green function would be 
confined to narrow regions at the left and right edges of the plot. For the 
case a z = 1.2 m, it would be confined to extremely narrow regions at the left 
and right edges. 




-0.08 -0.06 -0.04 -0.02 0.02 0.04 0.06 0.0S 0.02 0.04 0.06 0.08 0.] 0.12 0.14 0.16 



Figure 1: Left: On-axis charge density and Green function vs. z, as they appear in physical 
space, for the case cr z — 1.2 cm. Right: On-axis charge density and Green function as 
stored in memory. 

Figures 2, 3, and 4 show convolution results for the three different bunch 
lengths, a z = 1.2 cm, a z = 12 cm, and o z = 1.2 m, respectively, for grid 
sizes 64 x 64 x 128 up to 512 x 512 x 1024. The left hand side of each figure 
shows plots of the potential as a function of z on-axis for various grid sizes, 
comparing results based on the ordinary Green function and the integrated 
Green function. The right hand side shows the relative error of the calculated 
potential. In Fig. [2j a z is less than the pipe transverse size (1.2 cm vs. 4 cm); 
both the ordinary Green function and the IGF are accurate to better than 
1% for all the grid sizes shown. In Fig. [3j a z is somewhat larger than the pipe 
transverse size (12 cm vs. 4 cm); when the grid is coarse, the ordinary Green 
function has significant errors (a few percent to several tens of percent), while 
the IGF accuracy is 1% or better. In Fig. |1J a z is much larger than the pipe 
transverse size (1.2 m vs. 4 cm); in this case when the grid is coarse the 
ordinary Green function results exhibit huge errors (more than 100%), while 
the IGF accuracy is still 1% or better. As mentioned above, the accuracy 
of the IGF results is controlled by how well the grid resolves just the charge 
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density. For the non-IGF results, the accuracy depends on resolving both 
the charge density and Green function, and, due to the exponential fall-off 
of the Green function, a coarse grid gives unusable results. The relative 
error in the potential is plotted on the right hand side of the figures. These 
were obtained be plotting (4>-(fihighres)/(phighres), where <j) highres is the highest 
resolution result, obtained using the IGF with a 512 x 512 x 1024. 
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Figure 2: Left: On-axis potential vs. z showing the ordinary Green function result and 
the Integrated Green function (IGF) result for various grid sizes. The bunch is a Gaussian 
distribution with a x = a y = 6mm, cr z = 1.2 cm. Right: Relative error of the on-axis 
potential vs. z for grid sizes 64 x 64 x 128, 128 x 128 x 256, and 256 x 256 x 512. 




Figure 3: Left: On-axis potential vs. z showing the ordinary Green function result and 
the IGF result for various grid sizes. The bunch is a Gaussian distribution with a x = a y = 
6mm, <7 Z = 12 cm. Right: Relative error of the on-axis potential vs. z. 

Fig. [5] shows the ordinary Green function and the IGF on-axis in the 
range z G [—8 cm, 8 cm] for the three cases a z = 1.2 cm, a z = 12 cm, and 
a z = 1.2 m, all computed using 256 longitudinal grid points. The left hand 
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Figure 4: Left: On-axis potential vs. z showing the ordinary Green function result and 
the IGF result for various grid sizes. The bunch is a Gaussian distribution with a x = a y = 
6mm, a z — 1.2 m. Right: Relative error of the on-axis potential vs. z. 



side shows the ordinary Green function. The exponential fall-off of the Green 
function is obvious. Note that the ordinary Green function is the same for all 
3 cases, it is just sampled differently. Note also that for a z — 1.2 m, the grid 
is so coarse that only the central point is significant, all the others are so far 
into the region of exponential fall-off that they are effectively zero. If p were 
constant inside of a cell (which is approximately true for a z = 1.2 m), use of 
Eq. ^ would correspond to approximating the area under the true Green 
function with the area under the triangular shaped region of Fig. [5j left; so 
it is no surprise that the non-IGF result in Fig. [4] is highly inaccurate, with 
errors exceeding 100%. The right hand side of Fig. [5] shows the IGF on-axis. 
As before, the exponential fall-off of the Green function is obvious. However, 
it does not matter that the coarse data poorly resolves the IGF, because the 
accuracy does not depend on that, it only depends on having a fine enough 
grid that the charge density is well approximated by a linear function within 
a cell. For a z = 1.2 m, this approximately true even with as few as 128 grid 
points, hence the IGF result in Fig. [4] has good accuracy. 

Lastly note that, for a z = 1.2m with 256 grid points, the IGF is approx- 
imately equal to zero for all but the central point. Hence, the four terms in 



Eq. (10) for the discrete convolution/correlation solution of the Poisson equa- 
tion, such as the term on the left-hand side of Eq. ([6]), could be approximated 
by keeping just k' = in the sum over k', 



^max 3max 



h x h y h z pi/j' t kGi-i',j+j',o- (17) 

i'=i j'=i 
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Figure 5: Left: Ordinary Green function vs. z, for three different bunch lengths, using 
256 grid points. Right: Integrated Green function vs. z. 



In other words, the problem has been converted to one involving multiple 
2D - instead of 3D - mixed convolutions and correlations. The would reduce 
execution time by a factor equal to log2 of the number of longitudinal grid 
points. More generally, depending on the grid size and the exponential fall- 
off, it would be sufficient compute, 

Cat fmax neighbor s 

h x h y h z 2_^2^ piiji^Gi-iij+j^k-k'i (18) 

i'=l j'=l k'=k 

where N neighbors = if the IGF, , is appreciable only when k = 0, or 

^neighbors = 1 if the IGF is appreciable at k = and k — ±1, etc. 

5. Discussion and Conclusion 

A new method has been presented for solving Poisson's equation in an 
open-ended rectangular pipe. Compared with the Hockney method for iso- 
lated systems (Eq. ([5])) which can be computed with a single FFT-based 
convolution, the new method involves 4 mixed convolutions and correlations 



(Eq. (10)), where the convolutions and correlations are mixed among the 
problem dimensions. Starting with the Green function for a charge in an 
open-ended rectangular pipe (Eq. ([9])), an Integrated Green function (IGF) 
was derived (Eq. (Jl6|)). Simulations of a Gaussian beam in an open-ended 
pipe showed that the IGF approach is much more robust, i.e., it retains 
much better accuracy, than the non-IGF method over a wide range of bunch 
lengths. This is because the accuracy of the IGF approach depends only on 
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having a fine enough grid resolve the spatial variation of the charge density. 
In contrast, the non-IGF approach is sensitive to disparities between the spa- 
tial variation of the Green function and the charge density. Since the effort 
to compute the IGF is not much more than for the ordinary Green function, 
this argues that one should always use the IGF. 

In theory the calculation of the IGF, which is represented as an infinite 



series in Eq. ( 16 ), could make the simulation much more time consuming than 
the case for isolated boundary conditions. But in practice this is unlikely 
since simulations with isolated boundary conditions often re-grid at every 
time step (or as needed) to take account of the changing beam size. This 
would not be the case for the rectangular pipe Poisson solver if the beam 
filled most of the pipe transversely, since the IGF would be computed once, 
Fourier transformed 4 ways, stored, and reused. A possible exception would 
be if the longitudinal size were changing rapidly with time, but that is not 
usually the case. In fact, in many applications involving circular accelerators 
the longitudinal bunch oscillations are very slow. It should also be pointed 
out that it is not always necessary to compute the Green function over the full 
transverse cross section. If the beam happened to be very small transversely 
(but long compared with the transverse pipe dimensions so that the open 
treatment would be inappropriate), then the transverse physical mesh could 
be a small area centered on the z-axis, which would be much more efficient 
and accurate than gridding the full cross section. 

To implement this approach one needs to be able to control the direction 
of the FFTs in all dimensions. Unfortunately most packages don't offer this 
capability for built-in multidimensional FFTs (an exception was the Connec- 
tion Machine Scientific Software Library), so it is up to the programmer to 
implement it using ID FFTs. Such a capability has other uses. For example, 
when dealing with the Vlasov equation, 

^ + (?-fy)/-(V*-ty)/ = 0, (19) 

a split-operator, FFT-based method for advancing the distribution in phase 
space one time step can be expressed as 



/(<?, p, t) = e -& d *>e t< r*- d *>e-& d df(q, p, 0). (20) 

This can be evaluated easily by alternately performing forward and backward 
FFT's separately in q and p. Applications such as these argue that FFT 
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packages with built-in multi-dimensional capabilities should allow the user 
to control the FFT direction in every dimension separately. 
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Appendix A. FFT-based Convolutions and Correlations 

Discrete convolutions arise in solving the Poisson equation, as well as 
in signal processing. In regard to the Poisson equation, one is typically 
interested in the following, 



where x corresponds to the free space Green function, r corresponds to the 
charge density, and p corresponds to the scalar potential. The sequence {pj} 
has J elements, {r&} has K elements, and {x m } has M = J+K — l elements. 

One can zero-pad the sequences to a length N > M and use FFTs to 
efficiently obtain the {pj} in the unpadded region. To see this, define a 
zero-padded charge density, p, 




K-l 



j = 0,...,J-l 
, k = 0,...,K- 1 
j-k = -(K-l),...,J-l 



(A.l) 




(A.2) 



Define a periodic Green function, £ m , as follows, 




.TO 



if m = —(K — 1), . . . , J — 1 
if m = J, . . . , N — K, 
for any integer i 



(A.3) 
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Now consider the sum 



N-l 



N-l 



N-l 



1 W-> k (Yl PnW nk )(J2 ^W mk ) 0<j<N-l, (A.4) 



fc=0 



n=0 



m=0 



where W = e 2m / N . This is just the FFT-based convolution of {pk} with 
{£ m }. Then, 



K-l N-l 



N-l 



(A.5) 



n=0 m=0 

Now use the relation 

N-l 



k=0 



W {m+n - j > k = Ndm+n-^N (i an integer). 



(A.6) 



k=0 



It follows that 



K-l 



y 3 



■n+iN 



<j < N- I. 



(A.7) 



n=0 



But ^ is periodic with period N. Hence 

K-l 

<Pj = 



r «0-« < 3 < N - 1. 



(A.8) 



n=0 



In the physical (unpadded) region, j e [0, J — 1], so the quantity j — n in 
Eq. (A.8) satisfies —(if — 1) < j — n < J — 1. In other words the values 
of £j- n are identical to Xj- n . Hence, in the physical region the FFT-based 
convolution, Eq. (A.4), matches the convolution in Eq. (A.l). 

As stated above, the zero-padded sequences need to have a length > M, 
where M is the number of elements in the Green function sequence {x m }. 
In particular, one can choose A^ = M, in which case the Green function 
sequence is not padded at all, and only the charge density sequence, {r^}, is 
zero-padded, with k — 0, . . . , K — 1 corresponding to the physical region and 
k — K, . . . , M — 1 corresponding to the zero-padded region. 

The above FFT-based approach - zero-padding the charge density array, 
and circular-shifting the Green function in accordance with Eq. (A. 3) - will 
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work in general. In addition, if the Green function is a symmetric function 
of its arguments, the value at the end of the Green function array (at grid 
point J — 1) can be dropped, since it will be recovered implicitly through the 



symmetry of Eq. (A. 3). In that case the approach is identical to the Hockney 
method [DUE]. 



Lastly, note that the above proof that the convolution, Eq. (A.4), is 
identical to Eq. (A.l) in the unpadded region, works even when W~i k and 
W mk are replaced by and W~ mk , respectively, in Eq. (A.4). In other 
words, the FFT-based approach can be used to compute 



K-l 



Pj 



^ ] r k x j+k 



k=0 



j = 0,...,J-l 
k = 0, . . . , K - 1 
j — k — —{K — 1), 



(A.9) 



., J-l 



simply by changing the direction of the Fourier transform of the Green func- 
tion and changing the direction of the final Fourier transform. 
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